Implementing the Keldysh formalism into the ab initio Gaussian Embedded Cluster 
Method for the calculation of quantum transport 
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We discuss the key steps that have to be followed to calculate quantum transport out of equilibrium 
by means of the ab initio Gaussian Embedded Cluster Method recently developed by the authors. 
Our main aim is to emphasize through several examples that, if a sufficiently large portion of the 
electrodes is included in the ab initio calculation, which does also incorporate an electrochemical 
potential difference — tJ-R = eV , there is no need to impose an electrostatic potential V drop 
accross the system. 
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The effort devoted to investigate electronic transport 
through molecular bridges and metallic nanocontacts has 
sharply increased in recent yearsEIBEl. In particular, a 
variety of ah initio methods are being deve 
ing to catch the essentials of these systemsE 
These methods share most of the ingredients [Landauer- 
Keldysh formalism for transport and density functional 
(DP) theory], but differ in their numerical implementa- 
tion, which translates into their computational efficiency 
and prcdictibility power. A critical problem not yet fully 
solved is the way these piethods actually deal with out- 
of-equilibrium transporttJ'EJ'Ej. Most of the approaches 
that are being proposed impose from the outset an ex- 
ternal and uniform electuistatic potential drop across the 
molecule or nanocontactQu'EI. The latter is justified when 
considering planar electrodes, which in most cases do not 
correspond to experimental electrode geometries. This 
choice simplifies the numerical implementation, but can 
condition the outcome of the calculation. In this work we 
discuss how the Keldysh formalism can be implemented 
into ab initio transport schemes and, in particular, into 
the Gaussian Embedded Clustp-pMethod (GECM) re- 
cently developed by the authorsliSEJ. Conclusive numer- 
ical evidence is presented which shows that, if a suffi- 
ciently large portion of the electrodes is incorporated 
into the ah initio calculations, there is no need to add 
an external electrostatic potential to the self-consistent 
potential. This provides a natural, consistent, way of 
solving the electrostatic problem for generic electrode ge- 
ometries. 

In a single-particle description, inherent to DP calcu- 
lations, one obtains the current / through the molecule 
or nanocontact at non-zero temperature and finite bias 
voltage V by integrating the transmission probability: 



/=-/ T{E,V)[f{E^tJiL)~ f{E + ^,R)]dE, (1) 



where f{E — ^l,r) is the Permi distribution for the left 
and right electrodes whose electrochemical potentials are 
denoted by /iL and /iR , respectively, and A*l — A'r = eV . 
The transmission coefficient T{E, V) depends both on 
temperature and the applied bias, and is given by 



T{E,V) = TrpLG^+^PRG^- 



(2) 



where hats denote operators. The Green function oper- 
ators and the gamma operators depend on energy and 
the applied potential (both arguments dropped for con- 



venience 
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The latter are defined as Pl,r — r ^ R' 
where Sl,r denotes the self-energy operator for the part 
of the right and left semi-infinite electrodes which is not 
included in the ah initio calculation. The Green function 
operator is in its turn defined as 



{E ± i5)i -F- tf^ - S^^^ 



G(±) = /, (3) 



where / is the unity operator. On the other handL the 
density matrix out of equilibrium is obtained fromE3 



J-oo L J 



(4) 



where G^g are the matrix elements in a non-orthogonal 
basis of the lesser Green function operator given by 



G' 



i&^ 



/(S-AiL)fL + /(S-Mfl)f J G(-). (5) 



Most of the technical difficulties implicit in the evaluation 
of thei-abawe expressions have been discussed in previous 
worksOutl'Ej. There is, however, an important conceptual 
issue that, in our opinion, remains unsolved and worth 
clarifying. In the presence of a voltage difference be- 
tween electrodes one might be tempted to impose from 
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FIG. 1: Average onsite energies of the atomic orbitals versus 
atomic position in M19 — M4 — M19 nanocontacts (two 19 
atom (lll)-planes plus a four atom chain, see inset). The 
results correspond to V=Q (empty circles) and 5 V (filled 
circles) and M =A1 (a) and Au (b). 



the outset an external electrostatic potential drop across 
the molecule or nanocontact. As the potential profile de- 
pends strongly on the electrode geometry, one would have 
to solve Poisson equation with the boundary conditions 
appropriate for such geometry. A way around this wri 
lem is to assume a simple form for the electrodea3'l 
However, adding an external potential through the ho- 
mogeneous solution of the Poisson equation is unneces- 
sary if a significant part of the metallic electrodes has 
already been included in the initial cluster and an elec- 
trochemical potential difference is maintained. The dif- 
ference between left and right electrochemical potentials 
charges one electrode and discharges the other one, ef- 
fectively creating-, an electrostatic potential drop across 
the constrictiontZl. In general, once self-consistency is 
attained, an electrostatic potential difference V between 
atoms one or two layers inside opposite electrodes de- 
velops while they remain neutral. On the other hand, 
the potential difference between atoms on the surface of 
opposite electrodes is smaller than V since they carry 
the charges. In summary, we assume all the charges re- 
sponsible for the electrostatic potential across the con- 
striction to be part of our cluster instead of imposing 
a usually unknown external electrostatic potential. For 
a large enough number of atoms in the electrodes this 
should be essentially correct. 

The way we have actually implemented this self- 
contained, out-of-equlibrium scheme into the GECM ini- 
tially proposed in Ref. ^ and later fully developed in 
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FIG. 2: Transmission versus energy for V=0 (continuous line) 
and 5 V (dashed line) in the M19 - M4 - M19 nanocontacts 
of Figure 1, with M =A1 (a) and Au (b), are shown. 



Ref. is hereafter described. A first step consists in 
carrying out ra, DP calculation by means of the GAUS- 
S1AN98 codal3 of the region that contains the molecule 
or the set of atoms that form the contact between the 
electrodes plus a significant part of the electrodes. Once 
a reasonable accuracy has been achieved, the Green func- 
tion for the infinite system is calculated and the density 
matrix given by Eq. ^ used in the subsequent DP pro- 
cess. The infinite electrode selfenergies are calculated by 
means of the Bethe lattice approximation and a minimal 
basis set with pseudopotentials was used in the calcu- 
lations [see Ref. 15 for details]. In previous works, the 
DP calculations were carried out by means of the Becke 
three-parameter hybrid functional plus the|-eorrelation 
functional of Lee, Yang and Parr (B3LYP)N. A difii- 
culty specific to the out-of-equilibrium case comes from 
the fact that the density matrix in Eq. ^ is complex. This 
can be handled in two ways: i) treating exchange within 
DP, instead of carrying out a full Hartree-Pock calcula- 
tion (this can be easily done by using BLYP instead of 
B3LYP), and ii) incorporating the complex density ma- 
trix into the GAUSS1AN98 code. In the present case we 
have chosen the former as the second procedure jAmuld 
have required major changes in the standard codec2l. In 
order to get reasonably reliable results for finite bias, a 
root mean square error in the density matrix better than 
10~^ was mandatory. 

Calculations have been carried out for gold and alu- 
minum nanocontacts similar to those of the insets in Fig- 
ures 1 and 3. Two types of fee clusters have been con- 
sidered: (a) two (111) planes containing 19 atoms each 
plus a four atom chain labeled M19 — M4 — M19 (a ca- 
pacitor was described by taking out the four atom chain 
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FIG. 3: Same as Figure 1 for a (100) pyramidal arrangement 
(two pyramids containing 14 atoms each plus a four-atom 
chain; see inset) referred to as M14 — M4 — M14 (see text). 



leaving a wide gap between the planes); (b) two (100)- 
oriented pyramidal clusters plus a 4 atom chain labeled 
M14 - AM - MIA [(lOO)-oriented pyramidal nanocon- 
tacts as M5 - M5, M14 - M14 and M30 - M30 have 
also been considered]. Interatomic bulk distances have 
been taken for the whole cluster (2.86 A for Al and 2.88 
A for Au). All results correspond to zero temperature 
and an external bias voltage in the range 0-5 V. The 
structural stability of the system models is an important 
issue when one is interested in the interpretation of ex- 
perimental results. This is not the case here. 

Figure 1 shows the average on-site energies {5d6s6p for 
Au and 3s3p orbitals for Al) on all atoms of the nanocon- 
tacts M19 - M4 - M19 for zero and 5 V bias. This mag- 
nitude reflects only the electrostatic potential on each 
atom, not telling us anything about the chemical poten- 
tial (charge) on them. The on-site energies in the planes 
are less dispersed in the case of gold than in aluminum. 
This is probably due to the simple ,s-character of the 
wavefunctions at the Fermi level that gold has (as op- 
posed to the sp character in aluminum). At 5 V bias, the 
major drop in the potential occurs at the chain/plane 
contacts. While the total potential drop in gold (alu- 
minum) is 4.44 eV (4.6 eV) the drop between the first 
and the last chain atoms is only 1.67 eV (1.43 eV). In 
the case of zero bias, the results arc symmetric with re- 
spect to the geometric center of symmetry, as expected. 
Instead, a similar symmetry is absent for 5 V. Namely, 
whereas the potential drop between the left electrode and 
the first atom in the chain is 1.92 eV (1.69 eV) for gold 
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FIG. 4: Transmission versus energy for V=0 (continuous line) 
and 5 V (dashed line) in the M14 - M4 — M14 nanocontacts 
of Figure 3, with M= Al (a) and Au (b). 



(aluminum) nanocontacts, it is only of 0.85 eV (1.48 eV) 
between the chain end and the right electrode. This fea- 
ture seems common to all previously reported results and 
reflects not only bulk band structure features, but also 
details of the nanocontact geometry (see below). The 
transmission for V — and 5 V is shown in Figure 2. 
There are noticeable differences both in gold and alu- 
minum. In the latter the peak structure is significantly 
changed, whereas in gold the gap below -0.5 cV is par- 
tially filled. It is worth noting that although the total 
current calculated by integrating in a window [-V/2, V/2] 
the transmission for V=0 does not differ much from that 
obtained with the full non-equilibrium approach (more 
important changes are expected in complex molecules), 
the differential conductance is significantly different. 

The average on-site energies on all atoms of the 
nanocontacts M14 — M4 — M14 for zero and 5 V bias, 
are depicted in Figure 3. As above, the on-site energies 
in the outer planes are less dispersed in the case of gold 
than in aluminum. Now the potential drop along the 
chain is significantly different for the two metals: while 
in the case of gold the total potential drop is 4.29 eV 
and the drop between the first and the last chain atoms 
is only 0.162 eV, in aluminum these numbers change to 
4.56 and 2.043, respectively. These results support our 
claim in the sense that what matters are not only the 
bulk features of the metals but also the nanocontact ge- 
ometry. The asymmetry in the potential drop remarked 
above is also noted. Again, the differences in the trans- 
mission between zero and 5 V bias arc significant. The 
peak structure in aluminum is noticeably reduced and 
the gap in gold filled over a range larger than 1 eV. 

In Table I we report the fittings of the numerical re- 
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TABLE I: Fittings of the numerical results for the varia- 
tion with applied electrostatic potential V of the difference 
between the average on-site energies of the atomic orbitals 
at the outer planes in clusters similar to those of Figures 1 
and 2 (see text) and of current through them, with functions 
AE = aV'' and I/Go = cV^, respectively {AE, I/Go and V 
given in eV, Go being the quantum of conductance 2e^ /h). 



cluster 


a b 


c d 


A15 — A15 
AlU ~ AlU 

Am - Am 

AlU - AM - AlU 
A119 - AM - A119 
yl/19-vacumm-yl/19 


0.820 0.979 
0.882 0.992 
0.942 0.984 
0.900 1.008 
0.919 1.001 
0.953 1.001 


0.810 1.039 
0.807 1.08 
0.635 1.20 
1.196 0.865 
1.315 0.925 


Au5 - Au5 
AuU - AuU 
AuU " AuA - AuU 
Aul9 - Au4 - Aul9 
Aul9-Ya.cumm-Aul9 


0.803 0.935 
0.840 0.993 
0.882 0.971 
0.882 1.003 
0.928 1.001 


0.521 1.068 
0.668 1.067 
0.805 0.950 
0.886 0.91 



suits for the difference between the average energies of all 
atoms in the left and right planes in M19 - M4 - M19 
nanocontacts, and the two outer planes in the pyrami- 
dal nanocontacts, versus the applied external potential. 
The fittings were done with aV'', instead of assuming a 
linear relationship from the start. All regression coef- 
ficients were higher than 0.999. Although a and b are 
always rather close to 1, some significant features are 
worth of comment. There is a steady increase of a and b 
as the size of the outer planes increases. The values clos- 
est to one are found for the capacitors (the nanocontacts 
with a Ks 14A gap). In the latter case the coefficient 
a is only 7% (5%) smaller than 1 for gold (aluminum) 



nanocontacts, while the exponent cannot in practice be 
differentiated from one. The small deviations with re- 
spect to the " ideal" behavior are due to the finite charge 
that those planes carry and/or to a yet insufficiently large 
electrode. As expected, in pyramidal nanocontacts the 
results steadily approach the ideal behavior as the size of 
the cluster increases. These results constitute the main 
message of the present work. The results for the current 
versus voltage were also fitted with //Go = bV^ (see Ta- 
ble I). In this case the lowest regression coefficient is 0.99 
(in some cases the results are better fitted with second 
order polynomia). The significant deviations from a lin- 
ear behavior (Ohm's law) are not surprising in systems as 
small as those investigated here and have been reported 
previously. The current is in general higher in aluminum 
than in gold, which is consistent with the larger conduc- 
tance found for aluminum at the Fermi levelEj. Devia- 
tions from linearity are also more important in the case 
of aluminum, with the exponent increasing with the size 
of the pyramidal cluster. This is a consequence of a num- 
ber of facts, including the increasing contributions from 
the TT channel discussed in Ref. ^ 

Summarizing, we have discussed the practical proce- 
dures one should follow to implement the non-equilibrium 
Keldysh formalism into the Gaussian Embedded Clus- 
ter Method previously developed by the authors. In our 
view the most important outcome of this work is the 
conclusive evidence concerning the need of incorporat- 
ing a sufficiently large portion of the electrodes in the ah 
initio calculation. If this is done, the artificial addition 
of an applied external potential to the potential in the 
Schrodinger equation is unnecessary. 
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An error, recurrent in the literature, consists of writing the 
elements of the density matrix without the overlap matri- 
ces that appear in Eq. (5). This error is compensated by the 
erroneous writing of the Green functions, say for instance 



[{E + iS)S - Fno - 
where Slti are the matrix elements of the operator 



the retarded as G 
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E^^'' -I- S'^'' in a non-orthogonal (NO) basis. One can 
easily check by taking matrix elements in the operator 
equation (4) that, in a non- orthogonal basis, G 
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